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ABSTRACT 

We present new Keck images at 0.9 /im and OVRO 1.3 mm continuum im- 
ages of five Class I protostars in the Taurus star forming region. We analyze 
these data in conjunction with broadband spectral energy distributions and 8- 
13 iim. spectra from the literature using a Monte Carlo radiative transfer code. 
By fitting models for the circumstellar dust distributions simultaneously to the 
scattered light images, millimeter continuum data, and the SEDs, we attempt 
to distinguish between flared disks, infalling envelopes with outflow cavities, and 
combinations of disks and envelopes. For each of these circumstellar density 
distributions, we generate grids of models for varying geometries, dust masses, 
and accretion rates, and determine the best fits by minimizing the residuals be- 
tween model and data. Comparison of the residuals for best-fit disk, envelope, 
and disk+envelope models demonstrates that in general, models incorporating 
both massive envelopes and massive embedded disks fit the imaging+SED data 
best. The implied envelope infall rates for these disk-|-envelope models are gener- 
ally consistent with infall rates derived by previous investigators, although they 
are approximately an order of magnitude larger than inner disk accretion rates 
inferred from recent spectroscopic measurements. In addition, the disk masses 
inferred from our models are close to or larger than the limit for gravitationally 
stable disks, indicating that Class I disks may undergo periodic episodes of en- 
hanced accretion, perhaps as a result of gravitational instabilities. An important 
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caveat to these results is that in some cases, no single model can fit all of the 
imaging and SED data well, suggesting that further refinements to models of the 
circumstellar dust distributions around Class I sources are necessary. We discuss 
several potential improvements to the models, as well as new constraints that 
will become available with upcoming millimeter and infrared facilities. 

Subject headings: stars:formation — starsxircumstellar matter — stars:individual(lRAS 
04016+2610, IRAS 04108+2803B, IRAS 04239+2436, IRAS 04295+2251, IRAS 
04381+2540) — techniques: high angular resolution 

1. INTRODUCTION 

The canonical picture of low-mass star formation is that of a rotating, collapsing cloud 
of dust and gas that forms a protostar surrounded by a disk (e.g., Terebey et al. 1984; Shu 
et al. 1987, 1993). Different stages of this theoretical evolutionary process have been equated 
with observed differences in the spectral energy distributions (SEDs) of protostellar objects, 
which have been grouped into classes O-III based on their infrared spectral index and the ratio 
of sub-millimeter to bolometric luminosity (Lada & Wilking 1984; Lada 1987; Adams et al. 
1987; Andre et al. 1993). In this classification scheme. Class and I sources are still in the 
main accretion phase, and emit most of their radiation at far-IR and submm wavelengths due 
to reprocessing of light from the central protostars by dust grains in an infalling envelope. 
In contrast. Class II and III sources exhibit directly revealed pre-main-sequence stars in 
addition to emission from circumstellar disks. 

This classification scheme is an attempt to represent discretely a continuous evolutionary 
sequence, and there arc transition objects that create some blur between classes. Moreover, 
this sequence is defined solely from spatially unresolved SEDs, which contain only limited 
information about the circumstellar distributions. Thus, it is not clear that the observed 
differences in SEDs truly correspond to evolutionary changes in the circumstellar geometry. 
A crucial test is to constrain the geometry of material around members of these different 
evolutionary classes using spatially resolved images, and thereby either confirm or refute the 
evolutionary sequence inferred from spatially unresolved SEDs. 

Direct imaging of Class II sources has shown that the bulk of the circumstellar material 
lies in disks (e.g., McCaughrean & O'Dell 1996; Koerner & Sargent 1995; Eisner et al. 2003, 
2004), while there may be a small amount of material in tenuous envelopes (e.g., Grady 
et al. 1999; Semcnov ct al. 2004). Thus, Class II sources appear to be fully-assembled young 
stars surrounded by rotating disks from which they continue to accrete material. Direct 
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images of the less evolved Class I sources are relatively rare, due to the large extinctions to 
these embedded objects, and thus the circumstellar geometry for these objects is not well 
constrained. 

Studies of the emergent spectral energy distributions at wavelengths > 10 /im have 
provided important, albeit ambiguous, constraints on the circumstellar dust distributions 
for Class I objects. Models incorporating infalling, rotating, envelopes with mass accretion 
rates on the order of lO'^ Mq yr ^ are consistent with observed SEDs (e.g., Adams et al. 
1987; Kenyon et al. 1993a), and the compatibility of the derived mass accretion rates with 
statistically-inferred ages of Class I sources (Adams ct al. 1987; Myers ct al. 1987; Benson 
& Myers 1989; Kenyon et al. 1990) provides further support for these envelope models. 
However, for some objects whose SEDs can be explained by spherically-symmetric dust 
distributions, it has been suggested that nearly edge-on flared disk models may also be able 
to reproduce the observed SEDs (Chiang & Goldreich 1999). Observations that spatially 
resolve the circumstellar emission are necessary to remove the ambiguities inherent in SED- 
only modeling. 

The circumstellar geometry is constrained directly by spatially resolved images. Scat- 
tered hght at near-IR wavelengths (e.g., Kenyon et al. 1993b; Whitney et al. 1997, 2003b), 
with asymmetric emission morphologies observed in some sources, provides strong evidence 
that the circumstellar material is not spherically distributed around Class I objects. How- 
ever, since scattered light traces tenuous dust in the surface layers of circumstellar dust 
distributions, it does not allow an unambiguous determination of geometry or other circum- 
stellar properties. For example, while spherically-symmetric dust distributions can be ruled 
out, we show below that modeling of scattered light alone can not necessarily distinguish 
between a flattened envelope and a flared disk. 

Modeling of images at multiple wavelengths can provide tighter constraints on geometry, 
since emission at different wavelengths arises in different layers of the circumstellar material, 
with short-wavelength scattered light from low-density surface layers and longer-wavelength 
emission from deeper, cooler layers at larger radii. For the Class I source L1551 IRS 5, 
spatially resolved observations at several wavelengths (e.g., Strom et al. 1976; Keene & 
Masson 1990; Butncr ct al. 1991; Lay ct al. 1994; Ladd et al. 1995; Rodriguez ct al. 1998; 
Chandler & Richer 2000; Motte & Andre 2001) and detailed spectroscopy (White et al. 
2000) have been combined with SED modeling, providing important additional constraints 
on the models (Osorio et al. 2003). Moreover, the combination of low and high resolution 
millimeter observations facilitated distinction of compact disk emission and more extended 
envelope emission. Similar modeling of multi- wavelength observations for IRAS 04302-1-2247 
provided flrm constraints on the distribution of circumstellar material and the properties of 
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the circumstellar dust grains (Wolf et al. 2003). 

To further our understanding of the Class I population as a whole, we have obtained 
high spatial resolution observations of five additional Class I sources at multiple wavelengths 
which, when combined with SEDs, enable much tighter constraints on circumstellar dust 
models than available from any single dataset. We present images of scattered light at 0.9 
//m and thermal emission from dust at 1.3 mm. In addition, we augment SEDs compiled 
from the literature with new photometry at 0.9 fim, 18 fim, 1.3 mm, and 3 mm wavelengths. 
Using the three dimensional radiative transfer code MC3D (Wolf & Hcnning 2000; Wolf 
2003), we model our data in the context of three types of circumstellar dust distributions: 
1) rotating infalling envelopes; 2) flared disks; and 3) combinations of envelopes+disks. 

The best fits are obtained for models incorporating both an envelope and an embedded 
disk, although we show that pure disks or envelopes can reproduce certain aspects of our data. 
For each source, we discuss the properties of the best-fit models, including the geometry of 
the dust distributions and implied mass accretion rates. We use these results to help place 
Class I sources in the proper evolutionary context. Finally, we discuss refinements to the 
models that may improve agreement with observed multi-wavelength data, and describe new 
constraints that will become available from upcoming millimeter and infrared facilities. 



2. OBSERVATIONS 
2.1. The Sample 

Our sample consists of five Class I sources in the Taurus star forming region drawn from 
the larger sample of protostars studied by Kcnyon ct al. (1993a). This subset was selected 
based on three criteria: 1) sources must be detected at red optical wavelengths (-Fo.g^m > 10~® 
Jy), 2) this emission must appear resolved, and 3) sources must emit strongly at millimeter 
wavelengths (Fimm > 10^^ Jy). We initially observed the complete Kenyon et al. (1993a) 
sample at 0.9 iim with Keck/LRIS (see §2.3 and Appendix A), and chose only objects 
that satisfy our selection criteria for high angular resolution millimeter observations (§2.2) 
and detailed modehng. Our sample includes IRAS 04016+2610, IRAS 04108+2803B, IRAS 
04239+2436, IRAS 04295+2251, and IRAS 04381+2540. General properties of our sample 
are discussed here, and information about individual objects is given below in §4.1. 

The bolometric luminosities of our sample span 0.4 to 3.7 Lq, similar to the distribu- 
tion of luminosities for the complete sample of Kenyon et al. (1993a), which has a median 
luminosity of 0.7 Lq. The SEDs for our sample objects are similar in shape, exhibiting high 
extinction of the central objects, small amounts of mid-IR absorption, and peaks near 100 
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/im. While these SED characteristics are common to most other sources in the Kenyon et al. 
(1993a) sample, several of the SEDs for their sample exhibit deep 10 /xm absorption features. 
As discussed below, the depth of the 10 fim feature is correlated with source viewing an- 
gle; the dearth of objects with deep 10 /xm absorption in our sample may therefore indicate 
selection against edge-on sources. 

There arc other clear selection effects for our sample, since we are biased toward sources 
exhibiting bright scattered light and millimeter continuum emission. For example, by se- 
lecting sources exhibiting bright scattered light nebulae, we are selecting against face-on 
sources where the directly visible, point-like protostar would dominate the short-wavelength 
emission, as well as objects with extremely high line-of-sight extinctions, such as one might 
expect for massive spherical envelopes. Since edge-on sources produce dimmer scattered 
nebulosity, these may also be under-represented in our sample. Finally, our selection crite- 
ria for bright millimeter emission favors objects with higher masses of dust concentrated in 
smaller volumes. Future observations of larger samples with more sensitive instruments are 
necessary to investigate these potential biases. 



2.2. OVRO Observations 

We imaged each source in our sample at both 1.3 mm (230 GHz) and 3 mm (115 GHz) 
with the OVRO Millimeter Array between 2002 August and December. The 3 mm ob- 
servations were conducted in the "compact" array configuration, which provides baselines 
between 20 and 55 m, while the 1.3 mm observations were obtained mostly with the "equato- 
rial" configuration, providing baselines from 30 to 120 m. For IRAS 04295+2251 and IRAS 
04381+2540, 1.3 mm data were also obtained with the "high" configuration, which provides 
baselines from 35 to 240 m. Continuum data were recorded in four 1-GHz channels. We 
cahbrated the amplitudes and phases of the data using quasars near on the sky (< 20°) 
to our target sources. The flux of our target sources was cahbrated using observations of 
Neptune and Uranus, which have known millimeter fluxes, to calibrate the response of the 
instrument. All data calibrations were performed using the OVRO software package MM A 
(Scoville et al. 1993). Using the calibrated data for each source, we inverted the interfero- 
metric visibilities to create an image, and de-convolved and CLEANed the images using the 
MIRIAD package (Sault et al. 1995). We averaged the data using robust weighting (with 
a robustness parameter of 0.5) to obtain a good balance between sensitivity and angular 
resolution. 

From our data, we measured photometric fiuxes for our targets at both 3 mm and 
1 mm wavelengths. The emission is much stronger at 1 mm, as expected for thermally- 
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emitting dust, and thus the sources are more clearly detected at the shorter wavelength. In 
addition, since the angular resolution is inversely proportional to the observing wavelength, 
and because we obtained longer-baseline data at 1 mm than at 3 mm, our 1 mm images have 
better angular resolution. Our 1 mm images, shown in Figure 1, have angular resolution of 
~ 2", and rms sensitivity of ~ 1.5 — 2.5 mJy/beam. The measured fluxes, listed in Table 1, 
are integrated over the OVRO beam, which is ~ 2" at 1 mm, and ~ 6" at 3 mm. 

2.3. Keck/LRIS Observations 

We obtained images of our Class I sample at 0.9 /xm (Cousins /-band) using the W.M. 
Keck II telescope and the Low Resolution Imaging Spectrograph (LRIS; Oke ct al. 1995) in 
imaging mode. The observations were conducted on 1998 October 30-31 and 1999 December 
13. The integration time was 300 seconds, the seeing was 0''5 — 0''6, the field of view was 
6' X 8', and the plate scale was 0.21 arcseconds per pixel. Image processing included bias 
subtraction and flat fielding. The photometric fluxes of the target sources, integrated over 
a 6'.'3 (30 pixels) diameter aperture with background measured over an annulus extending 
from 26" to 30" (125-145 pixels), were calibrated using equatorial standards from Landolt 
(1992), assuming typical extinction coefficients for Mauna Kea (0.07 mag/airmass at / band; 
Beland et al. 1988); all sources were observed at an airmass of less than 1.1. 

All objects in our sample were detected at 0.9 /im, and show varying amounts of extended 
emission arising from scattered light due to dust-obscured central sources. Our images are 
shown in Figure 2. We determined absolute astrometry for these images by finding reference 
sources in common to our images and the 2MASS survey and computing a six coefficient 
fit for the image coordinates, plate scale, and distortion. The astrometry is accurate to 
~ 0"3, including the uncertainty in our astrometric solution and the overall uncertainty 
from 2MASS. 

Because we determine accurate positions for our /-band images, we can register the 
scattered light emission with the millimeter-wavelength emission observed with OVRO (§2.2). 
Since the longer wavelength radiation is optically thin, it traces dust mass and is therefore 
likely to be centered on the central source. In Figure 2, we have indicated the centroid of the 
millimeter continuum emission (and thus the likely position of the central protostar) with 
a cross. As we discuss below, the offset of the scattered light emission with respect to the 
central protostar provides valuable insights into source geometry and viewing angle. 
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2.4. Keck/LWS Observations 

Our sources were also observed at 10.3 and 17.9 fim with the Keck Long Wavelength 
Spectrometer (LWS; Jones & Puetter 1993) between 1999 August and 2000 December. LWS 
allows diffraction-limited (FWHM=0''2) imaging over a 10" field of view, although the ob- 
servations described here were obtained under poor seeing conditions, FWHM=0'.'3 — 0'.'6. 
With this seeing-hmited resolution, none of our targets are angularly resolved at 10 — 18 //m. 
Photometric fluxes were measured in a 1'.'92 (24 pixels) diameter aperture, with subtraction 
of the sky background measured in an annulus spanning 4''0-6''4 (50-80 pixels). Flat- fielding 
did not improve the data quality, and was therefore not applied before measuring photom- 
etry. Standard stars were observed each night to obtain atmospheric extinction curves, and 
curve-of-growth corrections were applied to each star to convert measured aperture photom- 
etry to infinite-aperture values. In addition, these Keck/LWS observations provided 8-13 
lim spectra, which have been reduced and analyzed by Kessler-Silacci et al. (2005). 

2.5. Spectral Energy Distributions 

SEDs were constructed using photometric fluxes from the literature and new measure- 
ments at 0.9 fim (§2.3), 18 fim (§2.4), 1.3 mm, and 3 mm (§2.2). Uncertainties are not 
available for much of the photometry from the literature, and therefore we adopt estimated, 
uniform photometric uncertainties in the modeling described below (§3). Since these uniform 
error bars lead to the over- weighting of higher fluxes, we consider the logarithms of the SEDs 
in our modeling. The fluxes used in our analysis, converted into units of Jy, are listed in 
Table 1 and plotted along with best-flt models below in §4. 

3. MODELING 

We use the three dimensional Monte Carlo radiative transfer code MC3D (Wolf & Ken- 
ning 2000; Wolf 2003) to model simultaneously the SEDs, 0.9 /xm scattered hght images, 
and millimeter continuum images of our sample of Class I sources (MC3D has been val- 
idated by benchmark comparison to other Monte Carlo and grid-based radiative transfer 
codes; Wolf et al. 1999; Pascucci et al. 2004). The basic model is that of a central emission 
source surrounded by obscuring dust of arbitrary geometry. The MC3D code solves the ra- 
diative transfer equation and determines the temperature distribution of circumstellar dust 
self-consistently. Moreover, the code takes into account absorption and multiple scattering 
events. The main inputs for this modeling are 1) parameters of the central star, 2) properties 
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of dust grains, and 3) geometry of circumstellar dust. We now describe which properties in 
these categories are varied in our modehng, as well as our assumptions for other parameters. 

Since our concern is primarily with understanding the geometry, we fix most of the 
stellar and dust parameters. For the models described below, we assume that the central 
object resembles a typical Class II T Tauri star. We therefore assume T* = 4000 K, — 2 
Kq, and M* = 0.5 M© (e.g., GuUbring ct al. 1998). Although few direct constraints exist 
on these parameters, spectroscopy of IRAS 04016+2610 has yielded an estimated effective 
temperature between ~ 3300 — 4200 K and a radius of ~ 0.9 R© (White & Hillenbrand 2004; 
Ishii et al. 2004), reasonably consistent with our assumptions. The dynamically-estimated 
stellar mass of IRAS 04381+2540 is ~ 0.2-0.4 (Brown & Chandler 1999), also compatible 
with our assumptions. 

Accretion onto the stellar surface may generate substantial luminosity in shocks near 
the stellar surface, which would supplement the stellar luminosity (Calvet & GuUbring 1998; 
GuUbring et al. 2000): ^central = L^, + Lace- In our modeling, we initially assume that 
-^central = -^*, but wc iucludc a scalc factor to allow for additional accretion luminosity (and/or 
slightly different stellar parameters from those assumed). As the last step in our modeling 
procedure, we find the value of ^central for which no scale factor is required. Although the 
spectral shape of emission from the accretion shock may differ from the stellar emission, we 
ignore this effect and assume that the temperature of all emission is 4000 K; the spectral 
shape is not critical given the extremely high optical depths in the inner regions of our 
models. In summary, all stellar parameters are constant from model to model, with the 
exception of a luminosity scaling that represents Lace- 

We assume that the dust grains are spherical with a power-law size distribution n{a) oc 
^-3.5 (^Mathis et al. 1977), where flmin = 0.005 fim and flmax = 1-0 yum, appropriate for ISM- 
like grains. We further assume that the dust is composed of a standard ISM mixture of 62.5% 
silicate and 25% ortho + 12.5% para graphite (e.g., Draine & Malhotra 1993; Weingartner 
& Draine 2001), with optical properties from Draine & Lee (1984). Previous authors have 
explored the effects of varying the chemical composition and particle size distribution of the 
circumstellar dust (e.g., D'Alessio et al. 1999, 2001). These parameters affect the overall 
SED, and 10 and 18 fim silicate features, to some extent. For simplicity, we keep the dust 
composition and particle size distribution fixed in our modeling. 

The main remaining input for radiative transfer is the density distribution, and we focus 
our efforts on exploring the parameter space associated with circumstellar geometry. We con- 
sider three classes of models: rotating infalling envelopes, flared disks, and envelopes+disks. 
Each of these models is defined from an inner radius, Rm, to an outer radius, -Rout- The 
inner radius is assumed to be 0.1 AU (comparable to the dust sublimation radius), although 
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this parameter is not crucial to the modeling given the high optical depths near to the cen- 
tral protostar. However, Rout may have substantial effects on models for the circumstellar 
material, and thus we allow this parameter to vary in our modeling. 

For each of the models considered in the following sections, several important parameters 
influence the density distributions, with concomitant effects on observed images and SEDs. 
In order to explore the effects of these parameters, we generate small grids of models by 
varying several important properties. Large grids are not possible due to the long nm-time 
of Monte Carlo radiative transfer codes (typically several hours per model), and we therefore 
consider only three values for each parameter. 

For the pure envelope model (§3.1), we vary the mass infall rate, M, the centrifugal 
radius, Rc, the outer radius, -Rout) and the inclination, i. The values of M are chosen to 
provide total envelope masses^ of 5 x 10^^, 10^^, and 5 x 10~^ Mq; since the envelope mass 
depends to some extent on parameters other than M (see Equation 3 and §3.1), the accretion 
rates explored include 27 discrete values between ~ 10~^ — 10~^ Mq yr~^. Sampled values 
of the other parameters are R^ = 30, 100, and 500 AU; R^nt = 500, 1000, and 2000 AU; 
and i — 5 — 90° in increments of 5°. For the pure disk model, we vary the disk mass, scale 
height, outer radius, and inclination; Mdisk = 10~^, 5 x 10~^, and 1.0 M©, ho — 5, 15, and 
25 AU, Rout = 500, 1000, and 2000 AU, and i = 5 - 90° in increments of 5° . Finally, for 
the envelope+disk model, we assume ho = 15 AU, but allow Mdisk to vary, in addition to 
varying several important envelope parameters: M, R^, Rout, and i. The range of envelope 
parameters are the same as for the pure envelope model, and we sample Mjisk = 10^'\ 10~^, 
and 1.0 Mq. All of the parameter values in these grids are chosen to bracket physically- 
plausible values. 

In order to compare with our observations, we use MC3D to compute self-consistently 
the temperature distribution, the SED from 0.1-5000 fim, the scattered light image at 0.9 /im, 
and the dust continuum image at 1.3 mm, for each model. The best-fit models are determined 
by minimizing the residuals between model and data for the combined SED-|-imaging dataset. 
Specifically, we minimize the quantity 

2_2_|_2_|_2 
Xr,tot — Xr,SED Xr,0.9/im Xr,lmm- K^) 

Here, Xr,SED reduced value of a model fitted only to the SED data, and Xr,o.9Mm 

and Xr,imm analogous quantities for the /-band and 1 mm images. Thus, x^.tot gives 

equal weight to the three datasets (as opposed to a true x? for the combined SED-|-imaging 



^All quoted envelope and disk masses increment the computed dust mass by an assumed gas to dust mass 
ratio of 100. 
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dataset, which would give more weight to the greater number of imaging datapoints relative 
to SED measurements). 

To calculate Xr values for each dataset, we estimate the uncertainties by assuming that 
the best-fit model (amongst all classes of models considered) for a given dataset has a reduced 
value of unity, and setting the error bars accordingly. This means that when a model 
is fitted to a given dataset (e.g., SEDs), the reduced is given by the sum of the squared 
residuals divided by the minimum value of this quantity for all of the models considered: 
denoting summed squared residuals by X"^, Xr = X^/min(X^). Thus, Equation 1 can be 
re-written as 

y2 

^-^'^^^ min{X^J ^ mm{Xl,^J ^ min(XL^) " ^ ' 

As mentioned above, when fitting the SED we minimize the difference between the logarithms 
of the model and data, to emphasize the shape of the SED over the peak value. In addition, 
we allow a scale factor between the model SED and the observed fluxes. When fitting 0.9 
yum and 1 mm images, we first rotate the observed images to a position angle where the 
brightest scattered light is south of the millimeter emission, consistent with the position 
angle definition in our models. 

Once the best-fit model out of these grids has been determined, we attempt to "zoom 
in" on the solution with more finely gridded values of inclination. During this zoom-in stage, 
we also vary Lccntrai to produce the correct normalization between the modeled and observed 
SEDs (without a scale factor). Because the central luminosity is somewhat degenerate with 
the total optical depth of the system (e.g., Kenyon et al. 1993a), minor adjustments to M 
and/or M^isk may be required when Lcentrai is varied. Thus, when zooming in on a solution, 
we vary i, I/centrai, and M and/or Mdisk simultaneously. 

The effects of various parameters on the SEDs have already been explored by previous 
investigators (e.g., Adams et al. 1987; Kenyon et al. 1993a; D'Alessio et al. 1999), and some 
work has studied the effects on scattered fight images (e.g., Whitney & Hartmann 1992; 
Whitney et al. 1997). We attempt to improve on this previous work by illustrating the 
effects of model parameters on our combined SED-|-imaging dataset for a broad range of 
circumstellar geometries, and discussing how degeneracies between fitted parameters can be 
broken. In the remainder of this section, we describe the rotating infalling envelope, flared 
disk, and disk+envelope geometries. In the following section, we present the best-flt model 
from each category for comparison to our observations. 
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3.1. Rotating, Infalling Envelope 

The density distribution for a rotating, infalling envelope is given by (e.g., Ulrich 1976; 
Cassen & Moosman 1981; Terebey et al. 1984), 

Penv(r, e) = ^{GMyy^/^ (i + J^) + 2^1^) . (3) 

Here, r is the radial coordinate, is the mass of the central star (assumed to be 0.5 M©), 
M is the mass infall rate, Rc is the centrifugal radius, fj, — cos 9 defines the angle above the 
midplane, and /Iq defines the initial streamline of the infalling material. 

For each value of r,9, there is a unique value of /^o, since the streamlines of infalling 
particles do not cross. Before determining the density distribution using Equation 3, we 
analytically solve for /iq : 

r cos sin^ 1 ^ /^^o /^\ 

Rc cos 6*0 — cos 6* 1 — jji/ jjiQ 

This is a cubic equation with three solutions. The correct solution is the one for which sin 9q 
has the same sign as sin^^, and for which the solution is real. In other words, a particle that 
starts in the northern or southern hemisphere never leaves that hemisphere. The analytic 
solution is thus ^ 

/^o = C - ^, (5) 

where ^ is defined as 

(1/3) 



{^IjirRl + V729/x2r2i?4 + lQmi{r - R, 



^ .3 



^~ 2(V3)3i?, ■ 

We include an outflow cavity in our envelope models by decreasing the density by some 
factor, /cav, in the polar regions of the circumstellar distributions. The cavity shape is defined 

by 

p(r, Z> Zo + r^)^ /cav X Penv{r, z). (7) 

Here, z is the height above the midplane, Zq describes how close to the star the outflow 
cavity begins, and C, describes the opening angle and shape of the outflow. We fix = 1 
AU, since this parameter does not have major affects on the observed SEDs or scattered 
hght images. We also assume that within the outflow cavity, the density decreases by a 
factor of /cav = 4 relative to the rest of the envelope; this value of /cav produces the correct 
amount of extinction to match the observed short-wavelength SEDs. The cavity profile has 
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some effect on tfie structure of modeled scattered light emission. However, these effects are 
small compared with the other envelope properties; thus, we fix the value of C at 1.1 in our 
analysis. Using other types of cavities (e.g., cones) may have small effects on the scattered 
light images, but will have negligible effects on SEDs (e.g., Whitney et al. 2003b). 

We fit this envelope model to our data by varying M, Rc, the outer radius of the 
envelope, /?out- and the viewing angle of the model, i. Some of these parameters produce 

degenerate effects on the models, which is why we generate grids where these parameters are 
varied simultaneously. However, we address the effects of each parameter below in order to 
shed light on how our different datasets constrain various envelope properties. 

The accretion rate, Af , is one of the most important parameters in terms of its effects on 
the SED, since it is linked to the total envelope mass (Equation 3), and thus governs directly 
the optical depth of the dust surrounding the protostar (Figure 3). For higher mass accretion 
rates, the optical depth will increase, and less of the short-wavelength emission will escape 
(either directly or as scattered light). Moreover, the mid-IR absorption deepens, because 
of increased extinction of the stellar emission and absorption in the 10 and 18 /im silicate 
features. The short-wavelength radiation is re-processed and emitted at longer wavelengths, 
leading to more emission at far-IR through millimeter wavelengths for higher values of M. 
The higher optical depths associated with large mass accretion rates also push the visible 
scattering surfaces outward, leading to somewhat larger /-band images. However, for very 
large accretion rates, the near-IR emission may be quenched altogether (Figure 3). 

The centrifugal radius, Rc, defines the radius in the model where material falling in 
from outer regions joins rotationally-supported material in the inner region. Outside of 
Rc, the envelope is essentially spherical (except for the outflow cavity), while the density 
distribution is significantly fiattened interior to the centrifugal radius. In addition, matter 
builds up at R^, leading to a region of enhanced density. For smaller values of i?ci a larger 
fraction of the envelope material is spherically distributed, and the optical depth is larger 
for most sight-hnes. Thus, models with small Rc show large near-IR extinctions and deep 
mid-IR absorption (largely degenerate with the effects of M; Figures 3 and 4). In contrast, 
for larger R^., more of the envelope's mass is relegated to the midplane, making it easier 
for emission to escape. Thus, scattered emission is stronger for larger values of Rc- This 
parameter also has a critical effect on the modeled millimeter images: since the density is 
higher for radii < Rc, the size of the millimeter emission will correlate with the centrifugal 
radius (Figure 4). Thus, the sizes of our millimeter images constrain Rf. directly. 

The outer radius also affects the model, mainly by changing the optical depth; larger 
values of i?out lead to more mass in the envelope, and thus more extinction of short wavelength 
light and more emission of long- wavelength radiation (Figure 5). Smaller values of i?out may 
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require larger values of M and/or smaller values of Rc in order to correctly model the SED. 
-Rout is constrained strongly from our scattered light and millimeter images, since larger outer 
radii produce larger images (Figure 5). For our sample, the observed sizes of scattered light 
constrain i?out to be > 1000 AU, while the fact that the millimeter images are fairly compact 
rules out outer radii larger than ~ 1500 — 2000 AU. 

Inclination is a crucial parameter for modeling both the SED and the scattered light 
images. As noted by previous authors (e.g., Kenyon et al. 1993a; Nakazato et al. 2003) 
and shown in Figure 6, the SED changes drastically depending on viewing angle. At larger 
inclinations, where the observer's line of sight passes through more of the dense midplane, 
the overall optical depth of the model increases, leading to enhanced extinction at short 
wavelengths and absorption at mid-IR wavelengths. In contrast, for smaller inclinations the 
optical depth is lower, and for face-on models, the line of sight may go directly down the 
outflow cavity to the central protostar. In addition to affecting the amount of scattering 
or direct stellar radiation, the inclination also affects the shape of the scattered emission: 
for edge-on distributions, the scattered light arises mainly from the edges of the outflow 
cavity, producing a symmetric structure with two lobes. As the inclination decreases, one 
lobe brightens with respect to the other, and moves closer to the central protostar as the 
line of sight pierces further into the cavity. Finally, for inclinations close to face-on, the 
protostar is visible directly, swamping any scattered emission that might be present. Since 
our millimeter images are only marginally spatially resolved, inclination is not important in 
modeling this emission. However, the offset between 1 mm and 0.9 /im emission is sensitive 
to inclination; the visible scattering surface moves farther from the mass-sensitive millimeter 
emission for larger inclinations (Figure 6). 



3.2. Flared Disk 

The density distribution for a flared disk is given by (Shakura & Sunyaev 1973), 



where 




Here, r is the radial distance from the star in the disk mid-plane, z is the vertical distance 
from the mid-plane, is the stellar radius, and h is the disk scale height. We assume fixed 
values for the exponents of radial and vertical density distributions, a and respectively. 
For a flared disk in hydrostatic equilibrium, one expects /3 = 58/45 (Chiang & Goldreich 
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1997), and we adopt this value. Using the relation that results from viscous accretion theory, 
a — 3(/3 — 1/2) (Shakura & Sunyaev 1973), we obtain a = 2.37. These values for a and j3 
are similar to those used in previous modeling of circumstellar disks (e.g., D'Alessio et al. 
1999; Wood et al. 2002; Wolf ct al. 2003). We fit the disk model to our data by varying 
Mdisk, ^0, -Rout, and the viewing angle of the model, i. 

The disk mass affects our modeled SEDs, scattered light images, and millimeter images. 
Regarding the SEDs, higher dust mass increases the optical depth, leading to higher obscu- 
ration of the central star, deeper 10 and 18 silicate absorption, and enhanced re-emission 
at long wavelengths. As illustrated by Figure 7, the disk mass also has substantial effects on 
the scattered hght: for sufficiently small disk masses (< 10~^ M©) the central protostar is 
visible for inclination < 70° and this stellar radiation swamps any scattered light. However, 
for higher disk masses, the flared disk surface becomes sufficiently optically-thick to obscure 
the protostar, allowing scattered emission to be observed. The morphologies of millimeter 
images are only slightly affected by disk mass: more massive disks have millimeter emission 
detectable out to larger radii, leading to slightly larger observed images. 

The effects of ho are similar to the effects of Mdisk in some respects; however, SED-|-imaging 
data allows us to disentangle the two. Larger values of ho produce more absorption at short 
wavelengths, and if the scale height becomes large enough to obscure the central protostar, 
scattered emission becomes visible (Figure 8). These effects arc similar to those produced 
by higher values of disk mass. However, the effects on the longer-wavelength emission dif- 
fer; higher values of Hq increase the far-lR flux, but have little effect on fluxes or image 
morphologies at millimeter wavelengths. 

Changing i?out docs not significantly affect the SEDs, as long as the disk mass is held 
fixed (Figure 9). This is because the optical depth of the model is not affected substantially 
by -Rout- However, there are small changes due to the larger surface area of optically-thick 
dust for larger values of -Rout- In addition, larger outer radii produce larger images in both 
scattered light and millimeter emission, since disk material is distributed to larger radii. 

Disk viewing angle is a crucial property in modeling the SEDs and scattered light 
images (e.g., D'Alessio et al. 1999; Whitney & Hartmann 1992). Similar to the dependence 
of envelope model SEDs on inclination, more edge-on disk models exhibit deeper absorption 
at mid-IR wavelengths, and higher extinction of the central star (Figure 10). Inchnation 
also substantially affects the observed scattered light: for moderate inchnations {i < 55°), 
the central star is visible and dominates the /-band emission, while for larger inclinations 
the protostar is obscured and an asymmetric scattered light structure is observed. For 
nearly edge-on orientations, a symmetric, structure is observed, corresponding to the top and 
bottom surfaces of the flared disk. As for envelope models, our marginally resolved millimeter 
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images can not constrain well the inclination, although the offset between scattered light and 
millimeter emission does probe the orientation of disk models. 

3.3. Envelope+Disk 

In addition to pure envelopes and pure disks, we consider a model incorporating both 
an envelope and an embedded disk. The main difference between this model and the pure 
envelope model is an enhanced density in the midplane, especially interior to the centrifugal 
radius. As a result of this enhanced midplane density, more long-wavelength emission is 
produced than for a pure envelope model. Since the amount of long-wavelength emission 
for a pure envelope model depends primarily on the mass infall rate, M, additional long- 
wavelength emission from the disk component may allow good fits to the data with lower 
inferred values of M. The effect of the disk component on the SED is illustrated in Figure 
11. 

In addition to larger fluxes, inclusion of a disk component also leads to more centrally 
concentrated modeled scattered light and millimeter images. With our marginally resolved 
millimeter images, we can not probe such density concentrations directly. However, as we 
discuss further below, comparison of our interferometrically measured fluxes, which trace 
compact emission, with lower-resolution measurements (e.g., Motte & Andre 2001; Young 
et al. 2003) can provide some constraint on the relative amounts of compact and large-scale 
material, providing an additional test of how well these models fit the data. 

The envelope-|-disk model is implemented by computing the density distributions for 
both a disk and an envelope, and then setting the density of the combined model to be the 
greater of the two individual densities for a given position, (r, ^). The temperature distri- 
bution, spectral energy distribution, and images are then calculated self-consistently for the 
combined model using MC3D. In our current implementation of the disk+envelope density 
distribution, the outer radius of the disk cannot be specified independently of the outer ra- 
dius of the envelope. Thus, our models may not truly represent a physical disk+envelope 
model, where one might expect the disk component to end at Re- These models should 
therefore be regarded as only qualitative indicators that both disk and envelope components 
are needed to match the data. 
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3.4. Disk+Extinction Model 

Finally, we consider another variation of a disk+envelope model, incorporating a flared 
disk density distribution plus foreground extinction. The extinction in this model mimics the 
obscuring nature of the envelope in the disk+envelope model described in §3.3, but employs 
a simpler geometry. In the following sections, we refer to this variation as a disk+extinction 
model. 

We consider values of Ay ranging from to 60 mag, which allows foreground extinctions 
much higher than expected from ambient material in the Taurus region (e.g., Kenyon & 
Hartmann 1995). Extinctions higher than possible from ambient cloud material {Ay ~ 5) 
probably arise in material associated with the collapsing cloud cores surrounding the sources. 
Sub-millimeter fluxes observed in large-beam SCUBA maps (Young et al. 2003) can be used 
to estimate column densities, and thus extinctions toward our sample of Class I objects. 
Since higher extinctions are inferred for more compact dust distributions, models requiring 
larger values of Ay imply smaller dust structures. Based on the SCUBA fluxes (Young 
et al. 2003), extinctions larger than ~ 5 mag are only possible if the material is distributed 
primarily on scales smaller than ~ 5000 AU, on the order of expected envelope sizes. 

In this disk+extinction model, the morphologies of scattered light and millimeter con- 
tinuum images are unaffected by the uniformly-distributed obscuring envelope material, and 
thus the model images resemble those predicted by pure disk models. Since the dependence of 
the modeled images on various disk parameters is the same as for pure disk models (§3.2), we 
refer to Figures 7-10. However, the foreground extinction substantially alters the synthetic 
SEDs by ehminating much of the short-wavelength stellar flux observed for pure disks. This 
model thus provides a way to produce disk-like images at the same time as heavily reddened 
SEDs. 



4. RESULTS 

For each of the density distributions described in §3, we determine the model providing 
the best-flt to our combined SED+imaging dataset. The properties of the best-flt models, 
as well as the reduced residuals between models and data, are listed in Table 2, and 
the models are plotted in Figures 12-16. Our results clearly indicate that pure disk models 
are not applicable for our sample of Class I sources. Rather, the data suggest models 
incorporating a massive envelope with an outflow cavity, although most sources arc flttcd 
best by models including both envelopes and disks. In this section, we discuss general results 
of our SED and image fltting, then describe results for individual sources in detail in §4.1. 
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While the observed SEDs for our sample show fairly shallow absorption at mid-IR 

wavelengths (Figures 12-16), edge-on disk models produce deep mid-IR absorption (Figure 
10; see also D'Alessio et al. 1999; Wood et al. 2002; Wolf et al. 2003), In addition, our 
asymmetric scattered light images rule out disk models close to edge-on, since such models 
would produce symmetric structures (Figure 10). However, the extended scattered light 
images indicate that disk models viewed close to face-on are not applicable either, since 
such models would be dominated by emission from the point-like central protostar. Since we 
observe spatially resolved, asymmetric scattered light toward our sample objects, we know 
that we are not directly observing the star, as would be the case for face-on or moderately 
inclined disk models. 

Although there may be a narrow range of inclinations for which pure disk models can 
fit both the SEDs and scattered hght images of some of our sources, substantial foreground 
extinction (^20 mag) is needed to make such a model consistent with the small optical/near- 
IR fluxes for our sample. As indicated by Table 2 and Figures 12-16, disk-|-extinction models 
provide far superior fits than pure disk models with Ay = 0. However, the extinction values 
for the best-fit disk+extinction models are all much larger than expected from ambient ma- 
terial in the Taurus cloud (e.g., Kenyon & Hartmann 1995), indicating a large concentration 
of mass on small scales toward these Class I sources. Thus, the disk-|-extinction models 
imply the existence of massive envelopes in addition to disk components. 

Pure envelope models also fail to fit the data well for most targets in our sample (Table 
2). As illustrated in Figures 12-16, pure envelope models typically over-predict the peak 
fiux in the SED. As discussed further below, pure envelopes also predict more extended 
millimeter emission than actually observed. Thus, pure envelope density distributions also 
do not seem suitable for explaining most features of the data for our sample. 

For most, if not all, of the objects in our sample, the best fitting models (i.e., those for 
which the residuals between model and data are minimized) incorporate both an envelope 
and a disk (described either by disk-|-envelope or disk-|-extinction models; Table 2). However, 
while SEDs and 1 mm images are fitted best by disk-|-envelope models, the scattered hght 
images for some sources seem to favor pure disk or pure envelope models. This may suggest 
that our disk-|-envelope model is too simplistic, and that a more complex implementation of 
disk-|-envelope models would allow all of the data to be fitted simultaneously. Alternately, 
the dust grain size distribution or composition may differ from our standard assumptions or 
vary radially or vertically within the model density distributions. One piece of evidence in 
support of the former hypothesis is the fact that the best-fit pure disk or envelope models, 
which usually fit the scattered light images well, generally have larger values of i?out than 
the best- fit disk-|-envelope models (Table 2) . This suggests a flattened distribution of small 
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dust grains that extends beyond the outer radius of our best-fit disk+envelope models. We 
discuss this issue further in §5.1. 

One way to test the relative contributions of disks and envelopes is to compare the com- 
pact millimeter emission observed in our OVRO observations with more extended emission 
seen in lower-resolution single-dish observations (Motte & Andre 2001). As illustrated in 
the case of IRAS 04016-1-2610 in Figure 17, pure disks, pure envelopes, and disks-|-envelopes 
produce varying amoimts of compact and extended emission, and thus comparison of high- 
and low-resolution millimeter data can help to distinguish between models. As described in 
§4.1, our interferometric fluxes typically (with one exception) differ by ~ 20 — 50% from the 
single-dish measurements of Motte & Andre (2001) (Table 3), indicating that the sources we 
are observing produce most of their emission from compact inner regions. Thus, disks appear 
to be an important, if not dominant, component of the circumstellar distributions around 
our sample of Class I objects. Combined with the arguments presented above, this provides 
further evidence that the circumstellar material resides in some combination of disks and 
envelopes. 



4.1. Results for Individual Sources 

4.1.1. IRAS 04016+2610 

IRAS 04016+2610 is an IRAS source (Beichman et al. 1986) lying at the western edge of 
the L1489 dark cloud (Benson & Myers 1989), and driving a molecular outflow (Myers et al. 
1988; Terebey et al. 1989; Moriarty-Schieven et al. 1992). Previous modeling of this source in 
terms of an infalling envelope model found M ~ 5 x 10~^ Mq yr~^, Rc ~ 50 AU, i ~ 45 — 65°, 
and = 3.72 Lq (Kenyon et al. 1993a; Whitney et al. 1997). An independent estimate 
of the inclination {i ~ 60°) was obtained from observations of a compact molecular outflow 
(Hogerheijde et al. 1998). Previous observations of scattered hght at near-IR wavelengths 
(Tamura et al. 1991; Whitney et al. 1997; Park & Kenyon 2002; Ishii et al. 2004) showed a 
similar morphology to our I band observations. However, our observations enable for the flrst 
time accurate registration of the scattered light to the position of the central source, which 
allows better constraints on the circumstellar dust distribution than possible in previous 
analysis. 

Our best-flt model for IRAS 04016+2610 incorporates both a rotating, infalling envelope 
and an embedded disk. This model flts all of our data well, including the 0.9 /xm image, 
1 mm image, and SED (Figure 12). The properties of this model are M = 6 x 10~^ Mq 
yr^^ i?c = 100 AU, i?out = 2000 AU, and Mask = 0.01 M©. The properties of the envelope 
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and disk components are consistent with previous analyses in the context of either model 
individually (Kenyon et al. 1993a; Whitney et al. 1997; Hogerheijde 2001). 

In order to correctly fit the total luminosity of the system, we require a large central 
luminosity, ^central ~ 4.7 Lq. This is substantially larger than the stellar luminosity inferred 
from near-IR spectroscopy (Ishii et al. 2004; White & Hillenbrand 2004), and suggests a 
large accretion luminosity. Since the total accretion rate determined from our fitting of 
disk+envelope models is not substantially higher than for other sources, we suggest that this 
additional luminosity may be generated in a more active accretion shock near the protostellar 
surface. 

For IRAS 04016+2610, the flux measured in our ~ 2" beam is approximately 40% of the 
value measured by Motte & Andre (2001) in an 11" beam. This is somewhat larger than the 
ratios predicted by our best-fit models, and may indicate that there is additional compact 
emission that we have not accounted for in the models. For example, if the disk is confined 
to radii smaller than Rc, as opposed to it!out as assumed in our models, more mass would 
be concentrated at smaller radii. Finally, we note that although the pure envelope model 
appears to fit our combined imaging+SED dataset almost as well as the disk+envelope model 
(Table 2), the ratio of compact to extended emission for the pure envelope model is much 
smaller than observed, providing further support for the disk+envelope model. 

4.1.2. IRAS 04108+2803B 

IRAS 04108+2803 is a 22"-separation binary system (e.g., Duchene et al. 2004) in the 
L1495 region, and IRAS 04108+2803B is the component that appears less environmentally 
evolved based on its spectral energy distribution; it emits the vast majority of the far-IR 
emission from the system. IRAS 04108+2803B also shows large scatter in photometric 
observations, indicating that it may be a variable star (Kenyon et al. 1993a). 

Previous modeling of this object in the context of infaUing envelope models found M ~ 
5 X 10-6 Mq yr-^ ~ 70 - 100 AU, i ~ 30 - 60°, and = 0.63 L©. (Kenyon et al. 1993a; 
Whitney et al. 1997). In contrast, Chiang & Goldreich (1999) showed that the SED of this 
object can be fit well by a fiared disk model with i = 65°, f3 = 1.2, and -Rout = 270 AU. The 
small fractional polarization (5.1%) relative to other Class I sources (~ 20%; Whitney et al. 
1997) was invoked as further evidence that the scattered emission from IRAS 04108+2803B 
arises in a disk rather than in the walls of an outflow cavity (Chiang & Goldreich 1999). 
Finally, IRAS 04108+2803B was modeled as a T Tauri star in a disk that is dynamically 
warped by a hypothetical stellar companion (Terquem & Bertout 1996). 
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Our results show that the data for this object cannot be fitted by a pure disk model, 
although the best fit is obtained for a model incorporating a massive disk in addition to 
an envelope (either disk+envelope of disk+extinction; Table 2; Figure 13). The best-fit 
disk+extinction model implies a disk with a mass of 0.6 Mq surrounded by an envelope 
providing 25 mags of foreground extinction. The best-fit disk-|-envclopc model for this source 
imphes M = 5 x 10"^ Mq yr-\ = 30 AU, i?out=500 AU, and i = 24°. The central 
luminosity is close to what one would expect for a T Tauri star, Lcentrai ~ 0.4 Lq. The disk 
mass of our best-fit model is 0.5 Mq (Table 2); however, as discussed in §5.2, this is likely 
an over-estimate. 

The ratio of compact to extended millimeter emission (0.7; Table 3) implies that a 
large fraction of the flux is generated by a compact disk component. Thus, this source 
appears to be extremely disk-dominated, supporting the hypothesis of Chiang & Goldreich 
(1999). On the other hand, a massive envelope component also appears necessary to fit the 
imaging+SED data, consistent with previous models (e.g., Kcnyon et al. 1993a) and recent 
Spitzer observations which attribute 15.2 /im CO2 ice absorption to a cold envelope (Watson 
et al. 2004). 

4.1.3. IRAS 04239+2436 

IRAS 04239+2436 is an IRAS source (Beichman et al. 1986) with a bright scattered 
light nebula (Kenyon ct al. 1993b), which has been tentatively associated with high- velocity 
molecular gas (Moriarty-Schieven et al. 1992). Previous modeling of SEDs and near-IR 
scattered light (separately) in the context of rotating, infalling envelopes, found M ~ 2 — 
4 X 10~^ Mq yr-\ i?c ~ 10 - 70 AU, i = 30 - 55°, and = 1.23 Lq (Kenyon et al. 1993a; 
Whitney et al. 1997). 

Although the SED for this source is fit best by a disk+envelope model, the combined 
dataset is fit comparably well by a disk+extinction model with Ay = 20 (Table 2; Figure 14). 
Thus, while there seem to be both disk and envelope components, the exact density profile 
of the source may not be matched exactly by either the disk+envelope or disk+extinction 
models. As we discuss further in §5.1, the fact that the disk+envelope model fails to represent 
accurately our 0.9 /im image may suggest that an additional flattened distribution of small 
dust grains extending out to ~ 2000 AU (not currently included in the model) is necessary to 
fit all of the data simultaneously. Assuming that our best-fit disk+envelope model provides 
an accurate tracer of the mass, the properties of this source are M = 3 x 10~^ Mq yr"-*-, 
Rc ~ 30 AU, i?out = 1000 AU, Mdisk = 0.5 Mq, Lcc„trai = 1-5 Lq, and i = 34°. However, as 
discussed in §5.2, this value of Mdisk may be substantially over-estimated. 
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In addition to our modeling, there are several arguments suggesting that IRAS 04239+2436 
is surrounded by both a massive envelope and an embedded disk. The ratio of compact to 
extended millimeter emission (Table 3) shows that a pure envelope model cannot fit the 
data, and thus there must be a compact disk component. On the other hand, recent Spitzer 
observations find a deep 15.2 fim CO2 ice absorption feature, which likely arises in a large, 
cold envelope (Watson et al. 2004). These arguments provide further support for a combined 
disk+envelope distribution of material around this source, and motivate future modeling that 
can fit all of the data simultaneously. 

4.14. IRAS 04295+2251 

IRAS 04295+2251 has appeared point-like in previous ncar-lR imaging observations 
(e.g.. Park & Kenyon 2002), and is tentatively associated with a molecular outflow (Moriarty- 
Schieven et al. 1992). Previous investigators analyzed SEDs and scattered light images 
(separately) in the context of infalling envelope models, and found M = 1 — 5 x 10~^ M0 
yr-^ i?c = 70 - 100 AU, i = - 30°, and = 0.44 Lq (Kenyon et al. 1993a; Whitney et al. 
1997). It has also been suggested that a pure disk model may be able to fit the SED for this 
source (Chiang & Goldreich 1999). 

Our modeling shows that the best-fit is obtained for a density distribution including 
both envelope and disk components (Table 2; Figure 15). Formally, the best-fit is obtained 
for a disk+extinction model including a 1.0 Mq disk surrounded by an envelope that provides 
20 magnitudes of extinction. The disk+envelope model also provides a good fit to the data, 
and the properties of this model are M = 4 x 10"*^ Mq yr"\ ~ 100 AU, i?out = 500 
AU, Mdisk = 1-0 Mq, /^central = 0.6 Lq, and i = 22°. As for other sources in our sample, 
the inferred values of McUsk for best-fit models are probably too high (§5.2). The ratio of 
compact to extended miUimeter emission observed for IRAS 04295+2251 (0.42; Table 3) is 
consistent with that predicted for the best-fit disk+envelope model, and implies a relatively 
disk-dominated density distribution. 

4.1.5. IRAS 04381+2540 

IRAS 04381+2540 is in the B14 region, and is associated with high-velocity molecular 
gas (Terebey et al. 1989). In addition, extended near-IR emission has been observed toward 
this object (Tamura et al. 1991). Previous modeling in the context of a rotating, infalling 
envelope found M = 5 - 13 x 10"^ Mq yr-^ i?e = 50 - 300 AU, i = 30 - 75°, and = 0.66 
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L0 (Kenyon et al. 1993a; Whitney et al. 1997). Moreover, the inchnation was estimated 
independently from observations of a compact molecular outflow to be ~ 40 — 70° (Chandler 
et al. 1996). 

Our modeling shows that the best fit is obtained for a density distribution incorporating 
both an envelope and a disk (Table 2; Figure 16). The properties of the best- fit disk+envelope 
model are M = 9 X IQ-^ Mq yr-\ = 30 AU, Mdisk = 1-0 Mq, i = 38°, and L^entrai = 
0.6 Lq. The disk+extinction model also provides a good fit, and may provide a more 
accurate representation of the 0.9 fim scattered light image. These results indicate that the 
circumstcllar dust distribution around IRAS 04381+2540 is dominated by a disk, although 
a massive envelope is also required to fit the data. However, the values of Mdisk in Table 2 
are probably over-estimated (§5.2) and the disk mass may be closer to ~ 0.2 M© for this 
object. 

A potential problem with the best-fit disk-|-envelope model is that it does not correctly 
reproduce the observed ratio of compact to extended millimeter flux (Table 3). Rather, 
this ratio indicates that the density distribution for this source is envelope-dominated. This 
may indicate that there is envelope material distributed uniformly over large spatial scales, 
leading to large amounts of extended flux being resolved out in our OVRO image. This 
is consistent with an extended structure extending to large radii from the central source 
(> 3000 AU) observed by Young et al. (2003), and may favor the disk+extinction model, 
which allows a large mass of dust to be distributed out to larger radii. Regardless of the 
exact density distribution, it seems that both disk and envelope components are necessary 
to fit the data for this source. 

The high derived accretion rate for IRAS 04381+2540, as well as the strong indication 
from Table 3 that the source is envelope-dominated, suggest that this object may be evolu- 
tionarily younger than other sources in our sample. Moreover, the mass accretion rate may 
be underestimated for this source: the mass of this object has been estimated dynamically 
from miUimeter spectral fine observations to be ~ 0.2 — 0.4 Mq (Brown & Chandler 1999), 
somewhat lower than our assumed value of 0.5 Mq, which implies that the value of M derived 
above may be underestimated by ~ 10 — 30% (Equation 3). 

5. DISCUSSION 

In §4, we fit our combined imaging+SED data for a sample of Class I sources with 
several models for the circumstcllar dust distribution: flared disks, collapsing envelopes, 
and combinations of disks and envelopes. Of the models considered, disk+envelope (and/or 
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disk+extinction) density distributions generally provide the smallest residuals between mod- 
els and data (Table 2). The properties of best-fit models are similar for the five sources in 
our sample, and for most parameters the spread in best-fit values is less than an order of 
magnitude. These tightly-clustered values are not surprising given the selection criteria for 
our sample and the similarity of the observational data for our targets (§2.1; Figures 12-16). 
However, these similar parameter values may also be due in part to finite parameter sampling 
and limitations in our models. 

In this section, wc discuss potential modifications to our disk+envelopc models that 
may be required to provide better fits to combined imaging+SED data, and examine the 
physical plausibility of derived model parameters. We also use our results to understand 
better the evolutionary stage of Class I sources, and to place them in context relative to the 
better-studied Class II objects. Finally, we discuss how new astronomical instruments will 
improve constraints on the circumstellar dust distributions for Class I objects. 

5.1. Lcirge-Scale Geometry 

While previous investigations of Class I sources argued for either pure envelope or pure 
disk models (e.g., Kenyon et al. 1993a; Whitney et al. 1997; Chiang & Goldreich 1999), our 
results based on simultaneous modeling of scattered hght images, thermal images, and SEDs 
(§4) indicate that the most suitable dust distributions likely incorporate massive envelopes 
and massive embedded disks. Given the large envelope centrifugal radii determined for our 
sample, the fact that massive disk components are also required is not surprising. In models 
of rotating infalling envelopes, the centrifugal radius demarcates the point at which infalling 
material piles up due to conservation of angular momentum (e.g., Terebey et al. 1984; Kenyon 
et al. 1993a). Because the centrifugal radius grows with time as the collapsing cloud rotates 
faster {Rc oc t^; e.g., Hartmann 1998), material should have previously piled up within the 
current value of Rc, creating a dense disk. Moreover, viscous spreading tends to smear out 
the material piled up at Rc into a more disk-like distribution. Since this material is not 
accounted for in the envelope density distribution (Equation 3), it is not surprising that the 
addition of a disk component improves the agreement between model and data. 

Although models incorporating both disks and envelopes generally yield the best fits to 
our combined SED-|-imaging dataset (Table 2), the exact density distribution is not firmly 
constrained since disk+envelope and disk-|-extinction models often provide fits of comparable 
quality. In addition, the scattered light images for some sources are fit better by pure disk 
or envelope models with large (~ 2000 AU) outer radii. Since scattered light probes trace 
material in surface layers, our 0.9 /^m images may contain contributions from flattened 
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distributions of tenuous material at larger radii. Moreover, our model assumes that the 
dust grain properties are the same everywhere. If one takes into account the effects of dust 
settling (e.g., Goldreich & Ward 1973), then the grains in the disk surface layer and outer 
envelope may be smaller than those in the dense midplane, allowing more scattering for a 
given mass of dust (e.g.. Wolf et al. 2003; Whitney et al. 2003a). 

The large polarizations observed toward Class I sources in near-IR polarimetric imag- 
ing observations support the notion of extended distributions of small dust grains (Whitney 
et al. 1997; Lucas & Roche 1998). At the edges of near-lR scattered light nebulosity ob- 
served for our sample, where most of the light is single-scattered at approximately 90°, linear 
polarizations are 70-80% (Whitney et al. 1997). In contrast, the maximum polarizations pre- 
dicted^ by our disk-|-envelope /-band model images are ~ 25%. Moreover, the integrated 
polarizations in our model images are substantially lower than the integrated near-IR polar- 
izations measured by Whitney et al. (1997), re-enforcing our hypothesis that there may be 
a population of small grains at large radii not included in our models. 

Since our best-fit models indicate that Class I sources are probably surrounded by 
envelopes with outflow cavities (in addition to disks) , we expect observed outflows from these 
objects to he in the middle of these cavities. For IRAS 04016+2610, IRAS 04239+2436, and 
IRAS 04381+2540, outflows have been observed at position angles of 165°, 45 — 60°, and 
10°, respectively (Hogcrheijde et al. 1998; Gomez et al. 1997; Saito et al. 2001). Comparison 
of these position angles with Figure 2 demonstrates that the outflows lie in the middle of 
the observed scattered light structures. For IRAS 04295+2251 and IRAS 04108+2803B, 
no geometrical information about molecular outflows is available. Thus, in all cases where 
outflow geometries can be derived, they are consistent with expectations from our modeling. 

The range of inclinations for best-fit models is small, with a spread of < 20° across all 
sources for a given model. This may be explained to some extent by selection effects: bright 
scattered light can only be observed if targets are inclined sufficiently to block the light from 
the central star. Moreover, scattered hght is somewhat brighter for moderate than for edge- 
on inclinations, which may bias our sample against edge-on sources as well. Nevertheless, 
the tight clustering of inclinations (Table 2) may be due in some part to limitations of 
our models in accurately representing the true density distributions. For example, for the 
disk+envelope model the assumed disk density profile may favor small inclinations for which 
the observer's line of sight passes through the fiared surface of the disk rather than the 
dense midplane, while inclinations > 30°, for which the line of sight does not pierce directly 



^The MC3D code has the capabihty of computing scattered hght images in different Stokes parameters, 
allowing recovery of polarimetric information (Wolf 2003). 
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down the outflow cavity, are required to reproduce the extended scattered hght structures 
and heavily reddened SEDs. Observations of larger samples and refinements to the models 
are necessary to provide a more reliable estimate of the inclination distribution of Class I 
sources. 



5.2. Disk and Envelope Masses 

As seen in Table 2, the masses of best- fit disk models, and the masses of the disk 
components of disk+envelope models, span a range of values from ~ 0.01 to 1.0 Mq. These 
high disk masses result naturally from the high optical depths of these models (r > 10^ — 10® 
at 100 fim, and > 10 — 1000 even at 1 mm). For pure disk models, large masses are required 
to extinct the central star (thus allowing scattered light to be observed) and to produce 
the observed millimeter emission (Figure 7). Moreover, the high optical depths that result 
from large disk masses lead to large amounts of cold dust in the midplane, which emits most 
of its radiation at millimeter wavelengths. For disk+envelope models, large disk masses 
enable enhanced millimeter emission without substantially altering the peak flux at shorter 
wavelengths. Moreover, substantial disk components lead to more centrally-concentrated 
millimeter emission, consistent with observations (Figure 17; Table 3). 

Disk masses larger than ~ 0.1 Mq may be unphysical, since they are gravitationally 
unstable (e.g., Laughlin & Bodenheimer 1994). Speciflcally, when the mass of a rotating disk 
is larger than approximately M^H/R (where H/R is the disk aspect ratio), the disk becomes 
gravitationally unstable and rapidly transfers angular momentum outward, resulting in rapid 
accretion of material onto the central protostar. This accretion process occurs on the order 
of the outer disk dynamical timescale (~ 10'^ — 10"^ yr), and is thus very fast compared to 
the inferred infall rate and should rapidly bring the disk mass down to a stable level. For 
our best-flt disk models, H/R ~ 1/10, and for a 0.5 M© star the stability criterion requires 
Miisk i$ 0.05 M0. While numerical simulations suggest that disks may remain stable with 
slightly higher masses, ~ 0.3M* Ri 0.15 M0 (e.g., Laughlin & Bodenheimer 1994; Yorke et al. 
1995), our models still require disk masses far in excess of this value. 

One implication of these un-physically-large disk masses is that rotating disk models 
seem untenable for modehng the imaging and SED data for Class I sources. Our massive, 
best-flt disk models may therefore resemble the magnetically-supported "pseudo-disks" pro- 
posed to occur during the early stages of protostellar collapse (Galli & Shu 1993a,b). Thus, 
while disk+cxtinction models provide good fits to the data for some sources (§4; Table 2), 
the large disks implied by these models may not correspond to the Keplerian disks seen at 
later evolutionary stages. 
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These large disk masses are also a concern for the best-fit disk+envelope models, where 
disk components are necessary to fit the SEDs well (Figures 12-16), and to provide the 
correct ratios of compact to extended millimeter emission (Table 3). As illustrated in Figure 
17, the most important contribution of the disk component occurs within the envelope's 
centrifugal radius, while outside of Rc, the millimeter emission from a disk-|-envelope model 
has a similar profile to the emission from a pure envelope. If the disk component is truncated 
at Rc (as opposed to i?out) as implemented in our radiative transfer code), the disk mass can 
be decreased by a factor of ~ 5, while the compact disk emission within Rc will still produce 
ratios of compact to extended millimeter emission consistent with observations. In order to 
produce the correct total millimeter flux, the envelope mass must be increased to compensate 
for the loss of emission from larger disk radii. As indicated by Figure 17, these outer disk 
regions contribute approximately 50% of the millimeter fiux at radii > Rc. By truncating 
the disk component of disk+envelope models at Rc, the disk masses listed in Table 2 may 
be decreased by a factor of ~ 5, while the envelope masses may increase by factors of ~ 2. 
Thus, true masses of the disk components for disk+envelope models may be < 0.2 M©, close 
to the limit of gravitational stability and therefore more physically plausible than the values 
hsted in Table 2. 

Even for disks truncated at Rc, the disk masses are often larger than the masses of the 
envelope components of disk+envelope models. Envelope masses range from 5 x 10~^ to 
0.05 Mq, and sources with smaller envelope masses tend to have relatively large disk masses, 
suggesting that different objects in our sample may be more or less disk-dominated. However, 
large-scale emission from extended material belies this trend to some extent: for example, 
while IRAS 04381+2540 appears to be the most disk-dominated source in our sample based 
on estimated disk and envelope masses, Table 3 indicates that there is probably a substantial 
extended dust component as well. Thus, it is hard to draw conclusions about the relative 
evolutionary stages of different objects in our sample based on inferred disk and envelope 
masses. 



5.3. Evolutionciry Stage 

An issue closely related to the properties of circumstellar material in Class I sources is 
their evolutionary stage. While fiared disk models are consistent with the SEDs of Class II 
objects (e.g., Kcnyon & Hartmann 1987; Chiang & Goldreich 1999; Dullcmond ct al. 2001; 
Eisner et al. 2004; Lcinert et al. 2004), we have shown that such models are not suitable 
for the Class I objects in our sample. Indeed, we have found that models incorporating 
both envelopes and disks provide a better match to the data. This finding would seem to 
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support the standard assumption that Class I sources are more embedded, still surrounded 

by massive envelopes, and potentially younger than Class II sources. On the other hand, 
spherically symmetric models (e.g., Larson 1969; Shu 1977) are not compatible with our 
data, indicating that Class I sources are at a stage intermediate to cloud cores and star+disk 
systems. 

However, other investigators have suggested that Class I and II sources actually may be 
at similar evolutionary stages. One piece of evidence in support of this hypothesis is that the 
spectroscopically-determined stellar ages of Class I and II sources appear indistinguishable 
(White & Hillenbrand 2004). It has also been suggested that differences in observed circum- 
stellar properties between the two classes are caused by different viewing angles, rather than 
different amounts or geometry of circumstellar material (Chiang & Goldreich 1999; White 
& Hillenbrand 2004). However, our results show that the five Class I objects in our sample 
are viewed at moderate inchnations (~ 30 — 50°), refuting earlier suggestions that optically- 
visible Class I objects arc predominantly viewed edge-on, and arguing that the circumstellar 
material around Class Is is indeed less evolved than that around Class lis. 

Another argument supporting similar evolutionary states of Class I and II objects is 
the fact that measured accretion rates, pertaining to the transfer of material from the inner 
disk onto the central star, appear similar for the two classes (White & Hillenbrand 2004). 
In contrast, the new data and modeling presented in this paper confirm earlier estimates of 
mass infall rates from the envelope onto the disk orders of magnitude higher than the derived 
inner disk accretion rates. The fact that disk and envelope accretion rates are not the same 
suggests periodic "FU-Ori" episodes where the accretion rate temporarily increases by more 
than an order of magnitude, likely due to a gravitational instability in the accretion disk 
(e.g.. Bell & Lin 1994; Hartmann & Kenyon 1996, and references therein). 

As noted by White & Hillenbrand (2004), the episodic accretion scenario would imply 
more massive disks in Class I sources relative to Class II objects. While Class I objects 
appear to be surrounded by larger total masses of dust, the emission on scales smaller than 
~ 12" indicates similar masses in the compact disk components (Andre & Montmerle 1994; 
Motte & Andre 2001; White & Hillenbrand 2004). However, the conversion of millimeter flux 
into mass depends on the optical depth, dust opacity and temperature; the higher optical 
depths and cooler dust temperatures in Class I sources may therefore lead to higher masses. 
Estimated disk masses for our best-fit disk+envelope models (§5.2) are larger than typically 
observed for Class II objects and are close to the values required for gravitational instability, 
providing support for the non-stationary accretion model. 
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5.4. Further Constraints on the Horizon 

We have demonstrated the power of combining spatially-resolved images at multiple 
wavelengths with broadband spectral energy distributions when modeling the dust distribu- 
tions around Class I sources. Low-density surface layers, hot inner regions, and cool outer 
regions are traced by optical/infrared scattered light, thermal mid-IR emission, and far-IR 
and millimeter emission, respectively. Because images at these various wavelengths probe 
different regions of the circumstellar material, and even different emission mechanisms, they 
place tight constraints on the range of circumstellar dust models consistent with the data, 
especially when used in conjunction. In this section we discuss the near-future prospects for 
spatially resolved imaging of Class I sources. 

Although our marginally-resolved millimeter images place important constraints on the 
centrifugal and outer radii (and viewing angle when combined with the position of near-IR 
scattered light), well-resolved images will provide better inclination constraints and useful 
additional information including direct measurement of the radial intensity profile and, in 
the best cases, the vertical intensity profile; taken together these constrain the density and 
temperature profiles. Analysis of spectral lines with different excitation conditions can probe 
directly the density and temperature structure. The enhanced angular resolutions and sensi- 
tivities of new and upcoming millimeter interferometers, including the SMA, CARMA, and 
ALMA, will enable detailed studies of Class I sources. 

At shorter wavelengths, the FWHM angular sizes of our typical disk-|-envelope model 
images (at 140 pc) are ~ 0'.'3 to 1''2 at 70 to 160 iim. Though smaller than the beam size 
of Spitzer, future far-IR interferometers such as the mission concepts SPIRIT and SPECS 
will resolve this thermal emission. At 10-20 /xm the emission is generated at small radii 
(< 10 AU), tracing dense material close to the star which can be studied with ground- 
based mid-IR interferometry. At 10 and 18 /xm, the FWHM sizes are ~ 30 and 40 mas, 
and the 1% emission contours lie at 140 and 280 mas, though at low fiux levels (< 30 
mJy). Current sparse-aperture interferometry with the Keck telescope can achieve angular 
resolutions of ~ 100 mas for sources brighter than a few Jy (e.g., Monnier et al. 2004), 
while future instruments like LBTI should yield additional improvements in resolution and 
sensitivity. Future capabilities such as mid-infrared instrumentation on TMT may also 
be able to achieve the required resolution and sensitivity. We note that different density 
distributions (e.g., pure envelopes or pure disks) can produce less centrally-peaked emission, 
which would lead to larger and thus more easily resolved images at infrared wavelengths. 

Finally, spectroscopy from the Spitzer Space Telescope will add valuable constraints on 
dust mineralogy and particle size distribution for Class I sources (e.g., Watson et al. 2004). 
Moreover, the shape of spectral features can provide a sensitive probe of dust optical depth 
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and source inclination. Since different molecules, ices, and dust species arise in different 
physical conditions, Spitzer observations will thus also provide powerful constrains on large- 
scale geometry. 

6. Conclusions 

We imaged a sample of five embedded Class I sources in the Taurus star forming region 
in 0.9 /im scattered light and thermal 1 mm continuum emission, and we analyzed these data 
together with spectral energy distributions and 10 //m spectra from the literature. Using the 
MC3D Monte Carlo radiative transfer code, we generated synthetic images and SEDs for four 
classes of models: 1) rotating infalling envelopes including outflow cavities; 2) flared disks; 
3) disks+cnvclopcs; and 4) disks+cxtinction. For each class of model, wc sampled small 
grids of relevant geometric parameters and determined the circumstellar dust distributions 
providing the best fits to our data in a sense. 

The imaging and SED data arc generally inconsistent with either pure disk or pure 
envelope models, and we find that the best fits are obtained with models incorporating both 
massive envelopes and massive embedded disks. Given the large centrifugal radii derived for 
our sample (i?c — 30—100 AU), the need for massive disks is not necessarily surprising. While 
not included in the rotating infalling envelope model, one expects a dense disk of material 
interior to Rc because of the growth of centrifugal radius with time: at earlier times, Rc was 
smaller, and thus material should have gradually piled up at successively larger radii out to 
the current value of Rc- Thus, models incorporating both envelopes and disks may be the 
most accurate (of those considered) for representing physical infalling envelopes. 

However, our results indicate that refinements to the models are necessary. For example, 
disk+envelope models including disk components truncated at Rc may yield disk masses 
closer to physically plausible values. In addition, our scattered light images point to the 
existence of extended, tenuous distributions of material, which must be included in the 
models in order to fit all of the data simultaneously. As more spatially-resolved data is 
collected over a broad range of wavelengths, further refinements to the models will likely be 
warranted. 

The overall geometry inferred for our sample of Class I objects is neither spherically- 
symmetric, as expected for the earliest stages of cloud collapse, nor completely fiattened as 
seen in Class II sources. Thus, our models confirm the picture where Class I sources are at 
an evolutionary stage intermediate to collapsing cores and fully assembled stars surrounded 
by disks. 
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The mass infall rates derived for our sample are between 3 — 9 x 10~^ Mq yr~^, consistent 
with previous results. However, these infall rates are more than an order of magnitude higher 
than accretion rates pertaining to the transfer of material from the disk onto the central star 
derived from high- dispersion spectroscopy (White & Hillenbrand 2004). This discrepancy 
argues for periodic "FU Ori" episodes of increased accretion: since infalling material is not 
accreted onto the star at the same rate, material piles up in the disk. Once the disk becomes 
unstable to gravitational instability, the accretion rate rises dramatically for a short period, 
depleting the disk mass and restoring stability. The "FU Ori" accretion hypothesis is also 
consistent with the high disk masses of disk+envelope models, which suggest disks near to 
the limit of gravitational stability. 

We have demonstrated the power of modeling broadband SEDs in conjunction with 
images at multiple wavelengths. In addition, we showed that 10 /im spectra can provide 
valuable additional constraints on circumstellar geometry, since the depth and shape of 
spectral features depends on dust optical depth and source inclination. Future observations 
with upcoming millimeter interferometers, the Spitzer Space Telescope, and other instruments 
will provide additional spatially- and spectrally-resolved information, greatly enhancing con- 
straints on the circumstellar dust around Class I objects. 
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A. LRIS Images of Larger Sample 

As discussed in §2.1, we obtained 0.9 /xm images using Keck/LRIS of the entire sample 
studied by Kenyon et al. (1993a). The details of these observations are described in §2.3. 
Although the analysis presented above concentrated on only those sources that satisfied our 
selection criteria (§2.1), we display the LRIS images for the larger sample in Figure 18 for 
completeness. 

Several of these objects are detected strongly at 0.9 //m: IRAS 04181+2654A, IRAS 
04248+2612, IRAS 04263+2436, IRAS 04264+2433AB, IRAS 04287+1802, IRAS 04303+2247, 
IRAS 04361+2547AB, IRAS 04368+2557, and IRAS 04489+3402. Four other sources ex- 
hibit faint smudges that may correspond to actual emission: IRAS 04169+2702, IRAS 
04181+2655, IRAS 04325+2402, and IRAS 04365+2635. The remaining sources were not 
detected, with a limiting /-band magnitude of ~ 24. Photometry for these objects, measured 
within a 6''3 diameter aperture, is listed in White & Hillenbrand (2004). 

Because these images are not very deep (< 300s), we can not comment on whether 
the pattern of detections and non-detections corresponds to variations in the underlying 
protostellar luminosities, different inclinations, or other factors. However, it is clear from 
Figure 18 that there are several sources exhibiting bright scattered hght structures; these 
sources are ideal candidates for future work. 



This preprint was prepared with the A AS IM^jX macros v5.2. 
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Fig. 1. — Images of our sample in thermal emission at 1.3 mm wavelength, obtained 
with OVRO. The origin of each image corresponds to the centroid of the millimeter emis- 
sion. The J2000 coordinates at the origin of the images are (04'^04'"43!11, +26°18'56:'5) 
for IRAS 04016+2610, (04^13"'54!72, +28°11'33^'0) for IRAS 04108+2803B, (04''26™56?29, 
+24°43'35:'l) for IRAS 04239+2436, (04'^32'^32?06, +22°57'26'.'3) for IRAS 04295+2251, and 
(04Ml°^12!71, +25°46'35^'3) for IRAS 04381+2540. 
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Fig. 2. — Images of our sample in scattered light emission at 0.9 /xm wavelength, obtained 
with Keck/LRIS. The images are registered with respect to the centroid of the millimeter 
emission for each source (Figure 1), indicated by a cross. 
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Fig. 3. — SEDs, scattered light images, and millimeter continuum images for a rotating, 
collapsing envelope with a range of mass accretion rates, spanning 10~^ to 10~^ Mq yr~^ (M 
increases from the bottom to top panels). This model assumes Rc — 30 AU, i?out — 1000 AU, 
and i — 45°. The /-band images were generated by convolving our model images with a 1" 
FWHM Gaussian, to simulate the seeing of our LRIS observations. Similarly, we convolved 
the 1 mm images with the actual beam that was used to observe IRAS 04016+2610. For 
both 0.9 /xm and 1 mm images, the contour intervals arc 20% of the peak flux, beginning at 
20%. The noisy appearance of some of the model images in this and subsequent flgures is 
due to the finite number of photons used in our Monte Carlo radiative transfer modeling, and 
does not reflect underlying substructure. As M increases, the strength of near-IR scattered 
emission decreases, although the shape of the emission broadens (note that for M — 10~^ Mq 
yr^^, the scattered emission is undetected in our models). Higher mass accretion rates also 
produce larger mid-IR absorptions. While the appearance of millimeter images is relatively 
unaffected by changes in M, the long- wavelength flux is enhanced for higher accretion rates. 
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Fig. 4. — SEDs, scattered light images, and millimeter continuum images for a rotating, 
collapsing envelope with a range of centrifugal radii, spanning 30 to 500 AU {Rc increases 
from the bottom to top panels). This model assumes M = 10^^ Mq yr~^, i?out = 1000 AU, 
and i = 45°. Larger values of -Rc lead to less spherically symmetric dust distributions, leading 
to increased absorption in the midplane, but less absorption out of the midplane. Thus, for 
large Rc, the central star becomes visible directly. The centrifugal radius also represents a 
region of enhanced density, and thus larger R^. lead to significantly larger millimeter emission 
regions. In fact for Rc — 500 AU, the millimeter emission is almost completely over-resolved 
in our OVRO observations. 
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Fig. 5. — SEDs, scattered light images, and millimeter continuum images for a rotating, 
collapsing envelope with a range of outer radii, spanning 500 to 2000 AU (-Rout increases 
from the bottom to top panels). This model assumes M — 3 x 10~^ Mq yr~^, — 30 AU, 
and i — 45°. As i?out increases, the amount of material in the envelope is increased (for a 
fixed M), increasing the optical depth and leading to shghtly higher absorption at near-to- 
mid-IR wavelengths. Not surprisingly, the images at both 0.9 fira and 1 mm become larger 
as i?out is increased. For the 1 mm images, the increased extended emission is manifested 
primarily by stronger sidelobe features, which are artifacts of the limited sampling in our 
intcrfcromctric observations (the 1 mm images have been convolved with the OVRO beam 
obtained for IRAS 04016+2610). 
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Fig. 6. — SEDs, scattered light images, and millimeter continuum images for a flared disk 
model at a range of viewing angles (i increases from the bottom to top panels). More edge-on 
models exhibit deeper absorption at mid-IR wavelengths, and higher extinction of the central 
star. For small inchnations {i < 30°), the central star is visible, and dominates the /-band 
emission. For moderate inclinations, an asymmetric scattered light structure is observed, 
while for nearly edge-on orientations, a symmetric, doublc-Iobcd structure is observed. The 
millimeter emission for this model is unresolved by our observations, and the model images 
are thus insensitive to inclination. 
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Fig. 7. — SEDs, scattered light images, and millimeter continuum images for a flared disk 
model with a range of disk masses (Mjisk increases from the bottom to top panels). For this 
model, i?out = 1000 AU, = 15 AU, and i = 60°. Higher disk mass increases the optical 
depth of the model, leading to higher absorption at short wavelengths and enhanced emission 
at long wavelengths. In addition, if there is sufficient dust mass, the central star is obscured, 
and scattered light is visible at short wavelengths. The flux in millimeter continuum images 
is increased for higher disk mass, and the emission becomes slightly more extended. 
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Fig. 8. — SEDs, scattered light images, and millimeter continuum images for a flared disk 
model with a range of disk scale heights {ho increases from the bottom to top panels). This 
model assumes M^isk = 5 x 10~^ Mq, Rout = 1000 AU, and i = 60°. Larger values of ho 
produce similar effects on the short-wavclcngth SED and scattered light images as larger 
values of Mdisk (Figure 7); higher absorption, and increased visibility of scattered light. 
However, the effects on the longer-wavelength emission differ; higher values of ho increase 
the far-IR flux, but have httle effect on miUimeter fluxes and lead to slightly more compact 
morphologies at millimeter wavelengths (corresponding to weaker sidelobe emission features) . 
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Fig. 9. — SEDs, scattered light images, and millimeter continuum images for a flared disk 
model with a range of outer radii (i?out increases from the bottom to top panels). This 
model assumes Mdisk = 5 x 10~^ M©, ho — 15 AU, and i — 60°. Since we have fixed the disk 
mass, different values of i?out do not significantly affect the opacity of the model. However, 
a larger outer radius produces a larger surface area of cool dust, which leads to a slight 
enhancement of the long-wavelength flux. Moreover, larger i?out will produce larger images 
in both scattered light and millimeter emission (the short-wavelength images shown here are 
dominated by light from the central protostar, and thus no scattered emission is visible). 
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Fig. 10. — SEDs, scattered light images, and millimeter continuum images for a flared disk 
model at a range of viewing angles {i increases from the bottom to top panels). This model 
assumes Mdisk = 5 x 10~^ Mq, ho ~ 15 AU, and i?out = 1000 AU. The effects of inchnation 
on disk models are similar to those shown for envelope models in Figure 6. More edge-on 
models exhibit deeper absorption at mid-IR wavelengths, and higher extinction of the central 
star. For moderate inclinations {i < 55°), the central star is visible, and dominates the I- 
band emission. For large inclinations, an asymmetric scattered light structure is observed, 
while for nearly edge-on orientations, a symmetric, double-lobed structure is observed. The 
millimeter emission becomes more extended in appearance as inclination increases, although 
this is a minor effect. 
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Fig. 11. — SEDs for envelope+disk models for a range of disk masses (Mdisk increases from 
the bottom to top panels). The inclination of the models has been fixed at 45°. The envelope 
component of the model has the following properties: M = 3 x 10~^ yr~^, M* = 0.5 
Mq, i?c = 100 AU, and i?out = 1000 AU. The various models illustrate the effect of adding 
progressively more massive (plotted bottom to top) disks to the model: Mdisk = 10~^, 10~^, 
and 1.0 M0. The main effect of the disk is on the SED long- ward of 100 /xm, where the 
dense, cool dust in the disk component adds substantial flux. The 0.9 /im and 1 mm images 
are also affected slightly, becoming more centrally concentrated for larger disk masses. 




Fig. 12. — SED data (Table 1), 10 fim spectrum (Kessler-Silacci et al. 2005), 0.9 /xm scattered 
light image, and 1.3 mm continuum image for IRAS 04016+2610, along with best-fit models 
for different assumed circumstellar geometries (Table 2). The observed images have been 
rotated on the sky by —10° (east of north) so that they have the same position angle definition 
as the models. Images have contour levels of 10% of the peak flux, beginning at 10% (~ 6(t) 
for the 0.9 fim images and 30% (~ 2.5a) for the 1 mm images. While pure disk models do not 
provide an adequate match to the data, extincted disks, pure envelopes, and disk+envelopcs 
match reasonably well. The disk+envelope, which provides the smallest residuals between 
model and data, is preferred. 
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Fig. 13. — SED data (Table 1), 10 /im spectrum (Kessler-Silacci et al. 2005), 0.9 /im scattered 
light image, and 1.3 mm continuum image for IRAS 04108+2803B, along with best- fit models 
for different assumed circumstellar geometries (Table 2). The observed images have been 
rotated on the sky by —90° (east of north) so that they have the same position angle definition 
as the models. Images have contour levels of 10% of the peak fiux, beginning at 30% 
(~ 2.50") for the 0.9 //m images and 30% (~ 3o") for the 1 mm images. Disk+envelope and 
disk+extinction models provide the best match to the data for this source. 




Fig. 14. — SED data (Table 1), 10 fim spectrum (Kessler-Silacci et al. 2005), 0.9 /xm scattered 
light image, and 1.3 mm continuum image for IRAS 04239+2436, along with best-fit models 
for different assumed circumstellar geometries (Table 2). The observed images have been 
rotated on the sky by —120° (east of north) so that they have the same position angle 
definition as the models. Images have contour levels of 10% of the peak flux, beginning at 
20% (~ 2.5a) for the 0.9 fim. images and 30% (~ 5a) for the 1 mm images. For this object, 
the best-fit to the complete data set is obtained for disk+cxtinction and disk+cnvclope 
models, although the SED and 10 /xm spectra are fitted best by the disk+envelope model. 
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Fig. 15. — SED data (Table 1), 10 \im. spectrum (Kessler-Silacci et al. 2005), 0.9 /xm scattered 
light image, and 1.3 mm continuum image for IRAS 04295+2251, along with best-fit models 
for different assumed circumstellar geometries (Table 2). The observed images have been 
rotated on the sky by 155° (east of north) so that they have the same position angle definition 
as the models. Images have contour levels of 10% of the peak flux, beginning at 10% 
(~ ISo") for the 0.9 ;um images and 30% (~ 5(t) for the 1 mm images. Disk+envelope and 
disk+extinction models provide a good match to the data for this source. Note that the 
modeled 1 mm images arc typically smaller than the observations in this case, since small 
values of R^. arc preferred by the SED and /-band data, forcing the best-fit model to a smaller 
centrifugal radius, and thus a smaller 1 mm image (see Figure 4). 




Fig. 16. — SED data (Table 1), 10 fim spectrum (Kessler-Silacci et al. 2005), 0.9 /xm scattered 
light image, and 1.3 mm continuum image for IRAS 04381+2540, along with best-fit models 
for different assumed circumstellar geometries (Table 2). The observed images have been 
rotated on the sky by 180° (east of north) so that they have the same position angle definition 
as the models. Images have contour levels of 10% of the peak flux, beginning at 30% (~ 2a) 
for the 0.9 /im images and 30% (~ 2.5a) for the 1 mm images. While pure disk models and 
pure envelope models do not provide an adequate match to the data, extincted disks and 
disk+cnvclopcs do. The disk+cnvclopc and disk+extinction models provide the smallest 
residuals between models and data (Table 2). 




Fig. 17. — Azimuthally- averaged radial profile of millimeter emission of best-fit models for 
IRAS 04016+2610 (Table 2). The pure envelope model (solid line) produces most of its 
millimeter emission near the centrifugal radius, 100 AU in this case. The pure disk model 
(dotted line) is more centrally-peaked, producing the majority of the millimeter emission at 
small radii. The disk+envelope model resembles the pure disk at small radii and the pure 
envelope at large radii. 



-54- 




Fig. 18. — Keck/LRIS images of the Class I source sample defined by Kenyon et al. (1993a). 
Each panel shows a 30" x 30" area. The LRIS images of the five sources analyzed in this 
paper are given in Figure 2, and those sources we did not analyze are shown here. Some 
of these objects are not detected in our observations, which may indicate different source 
luminosities or circumstellar dust geometries. 
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Table 1. Spectral energy distributions for our sample 



A (/um) 


04016+2610 


Source Fluxes, Fi, 
04108+2803 04239+2436 


(Jy) 

04295+2251 


04381+2540 


reference" 


0.55 


2.42e-05 










Kenyon et al. (1993a) 


0.63 


6.17e-05 










Kcnyon ct al. (1993a) 


0.79 


0.0001 










Kenyon et al. (1993a) 


0.89 


0.0001 


3.23e-06 


3.89e-05 


2.23e-05 


8.13e-06 


This work 


1.10 


0.0003 










Padgett ct al. (1999) 


1.22 


0.0055 


0.0004 


0.0010 


0.0055 


0.0005 


Kenyon ct al. (1993a) 


1.22 








0.0038 




Whitney et al. (1997) 


1.60 


0.0095 










Padgett et al. (1999) 


1.63 


0.0217 


0.0050 


0.0075 


0.0360 


0.0018 


Kenyon ct al. (1993a) 


1.63 


0.0233 


0.0242 


0.0076 


0.0229 


0.0029 


Park & Kenyon (2002) 


1.63 


0.0263 


0.0043 


0.0066 


0.0220 




Whitney et al. (1997) 


1.87 


0.0214 










Padgett et al. (1999) 


2.05 


0.0416 










Padgett ct al. (1999) 


2.19 


0.1160 


0.0226 


0.0421 


0.0943 


0.0101 


Kenyon et al. (1993a) 


2.19 


0.1045 


0.0420 


0.0464 


0.0690 


0.0147 


Park & Kenyon (2002) 


2.19 


0.1055 


0.0148 


0.0439 


0.0671 


0.0149 


Whitney et al. (1997) 


3.45 


0.4593 


0.0492 


0.2302 


0.1667 


0.0470 


Kenyon ct al. (1993a) 


4.75 


2.0465 


0.1702 








Kenyon et al. (1993a) 


7.80 


2.8 


0.31 








Myers et al. (1987) 


8.69 


2.10 


0.50 








Myers et al. (1987) 


9.50 


1.70 


0.49 








Myers et al. (1987) 


10.3 


2.50 


0.57 








Myers et al. (1987) 


10.6 


2.5682 


0.6020 




0.5883 


0.2087 


Kenyon et al. (1993a) 


10.7 


2.97 


0.76 


1.71 


0.72 


0.23 


Kessler-Silacci et al. (2005) 


11.6 


3.40 


0.44 








Myers et al. (1987) 


12.0 


3.8327 


0.7647 


1.7120 


0.6215 


0.4106 


Kenyon et al. (1993a) 


12.5 


4.90 


1.30 








Myers et al. (1987) 


17.9 


4.57 






1.07 




This work 


21.0 


10.6303 


2.6094 


4.4314 




0.1979 


Kenyon et al. (1993a) 


25.0 


16.3030 


3.6498 


6.9545 


1.8718 


2.7686 


Kenyon et al. (1993a) 


60.0 


48.1371 


7.6292 


15.2222 


3.4872 


10.5312 


Kenyon ct al. (1993a) 


100 


62.2772 


10.8225 


17.1525 


8.2097 


14.9392 


Kenyon ct al. (1993a) 


100 


55.5046 










Kenyon et al. (1993a) 


160 


46.6067 










Kenyon et al. (1993a) 


350 


4.7686 










Kenyon ct al. (1993a) 


450 


3.2926 


1.8094 








Kenyon et al. (1993a) 
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Table 1 — Continued 



A (/Km) 


04016+2610 


Source Fluxes, 
04108+2803 04239+2436 


(Jy) 

04295+2251 


04381+2540 


reference" 




9 81 QQ 








1 9000 


nugeiiieiju.e oc odiicieii yzuuvj 


450 


4.23 


1.13 


2.1400 


2.6600 


2.8199 


Young et al. (2003) 


800 


0.582 


0.085 


0.3330 


0.2410 


0.2890 


Moriarty-Schieven et al. (1994) 


800 


0.3143 


0.8859 


2.3303 






Kenyon ct al. (1993a) 


850 


0.5899 


0.1700 


0.4199 


0.4199 


0.5600 


Young et al. (2003) 


850 


0.43 








0.24 


Hogerheijde & Sandell (2000) 


1100 


0.180 


<0.100 


0.114 


0.094 


0.116 


Moriarty-Schieven et al. (1994) 


1100 


0.5317 


0.0882 


0.0768 






Kenyon et al. (1993a) 


1300 


0.150 


0.039 


0.170 


0.115 


0.300 


Motte & Andre (2001) 


1300 


0.0658 


0.0271 


0.0329 


0.0484 


0.0135 


This work 


2700 


<0.0068 








<0.0075 


Hogerheijde ct al. (1997) 


3000 


<0.0096 


<0.0210 


<0.0120 


<0.0210 




Ohashi et al. (1996) 


3000 




0.0058 


0.0120 


0.0087 




This work 


3369 


<0.0043 








<0.0050 


Hogerheijde ct al. (1997) 


3440 


<0.0180 




<0.0240 






Saito ct al. (2001) 


20000 


0.0005 










Lucas et al. (2000) 


35000 


0.0003 




0.0001 






Lucas et al. (2000) 


60000 


<0.0002 










Lucas et al. (2000) 



References. — "-Kenyon et al. (1993a) compiled SEDs from new data and the literature, and this work contains 
references for much of this photometry. 

Note. — Different observations were conducted with different angular resolutions, which explains some of the pho- 
tometric variations. In addition, all of the sources in our sample except IRAS 04295+2251 are known to be variable at 
short wavelengths (e.g.. Park & Kenyon 2002). 



Table 2. Best-fit Models 



Source 


Ar,tot 


Ar,SED 


v2 M 

Ar,0.9um Ar,lmm 


Rc 


Menv 


Mdisk 




Rout 


i 


-^central 










(M0 yr-i) 


(AU) 


(M©) 


(Mo) 


(AU) 


(AU) 


n 


(L©) 


Rotating Infalling Envelope Models 


IRAS 04016+2610 


4.39 


1.57 


1.40 


1.41 8 X IQ-^' 


100 


0.07 






2000 


36 


3.8 


IRAS 04108+2803B 


5.25 


2.90 


1.08 


1.26 6 X IQ-^ 


30 


0.05 






2000 


26 


0.7 


IRAS 04239+2436 


6.89 


3.93 


1.17 


1.77 8 X 10-6 


100 


0.07 






2000 


42 


2.3 


IRAS 04295+2251 


13.07 


10.05 


1.21 


1.79 9 X 10-6 


100 


0.01 






500 


34 


1.1 


IRAS 04381+2540 


4.25 


1.93 


1.01 


1.30 1 X 10-5 


100 


0.10 






2000 


30 


0.6 


Flared Disk Models 


IRAS 04016+2610 


14.30 


9.68 


3.35 


1.24 






0.10 


25 


2000 


52 


2.8 


IRAS 04108+2803B 


12.81 


10.30 


1.28 


1.21 






0.10 


25 


500 


52 


2.3 


IRAS 04239+2436 


11.21 


8.18 


1.73 


1.30 






1.00 


25 


500 


45 


1.0 


IRAS 04295+2251 


13.98 


10.60 


2.26 


1.11 






1.00 


15 


1000 


55 


0.9 


IRAS 04381+2540 


9.70 


7.56 


1.03 


1.10 






0.70 


15 


2000 


56 


0.5 


Disk+Envclopc Models 


IRAS 04016+2610 


4.12 


1.35 


1.51 


1.25 6 X 10-6 


100 


0.05 


0.01 


15 


2000 


37 


4.7 


IRAS 04108+2803B 


3.68 


1.40 


1.04 


1.22 5 X 10-6 


30 


5 X 10-3 


0.50 


15 


500 


24 


0.4 


IRAS 04239+2436 


4.73 


1.67 


1.72 


1.33 3 X 10-6 


30 


0.01 


0.50 


15 


1000 


34 


1.5 


IRAS 04295+2251 


5.83 


2.75 


1.93 


1.13 5 X 10-6 


30 


5 X 10-3 


1.00 


15 


500 


22 


0.6 


IRAS 04381+2540 


3.47 


1.34 


1.04 


1.08 9 X 10-6 


30 


0.01 


1.00 


15 


1000 


34 


0.6 


Disk+Extinction Models 


IRAS 04016+2610 


6.28 


1.82 


3.20 


1.25 




Ay = 10 


0.20 


25 


2000 


49 


5.3 


IRAS 04108+2803B 


3.34 


1.07 


1.05 


1.21 




Ay=25 


0.60 


25 


500 


17 


0.6 


IRAS 04239+2436 


5.01 


2.44 


1.21 


1.35 




Ay =20 


0.20 


20 


2000 


50 


5.4 


IRAS 04295+2251 


4.42 


1.08 


2.24 


1.09 




Ay =20 


1.00 


15 


500 


49 


1.2 


IRAS 04381+2540 


3.27 


1.12 


1.06 


1.07 




Ay =20 


0.70 


15 


500 


56 


1.0 



Note. — Best-fit models for different density distributions eonsidered in §3. Xi tot ^^'^ sum of the redueed chi-squared residuals between the 
data and model for the SED, 0.9 fim, and 1 mm images. The individual reduced chi squared values for each dataset are denoted by Xr SBD' 9;im' 
and Xr inim- t'^'^ ease of the disk+extinction models, entries in the Menv column correspond to foreground extinction. Although the disk masses 
listed in this Table are the results of our model-fitting, we believe that these are likely over-estimated (§5.2). 
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Table 3. Comparison of compact and large-scale millimeter emission 



Source 


-Rdata 


-Rdisk 


Renv 


-^disk+env 


IRAS 04016+2610 


0.43 


0.17 


0.06 


0.15 


IRAS 04108+2803B 


0.70 


0.60 


0.08 


0.57 


IRAS 04239+2436 


0.19 


0.53 


0.06 


0.29 


IRAS 04295+2251 


0.42 


0.28 


0.26 


0.35 


IRAS 04381+2540 


0.05 


0.12 


0.05 


0.32 



References. — In this table, i?data is the ratio of emission 
observed in our compact OVRO beam to the emission ob- 
served in the larger-beam survey of Motte & Andre (2001). 
^^disk, Rcnv, and i?disk+onv SLTC the ratios one would measure 
for pure disk, pure envelope, and disk+envelope models. For 
disk+extinction models (§3.4), the distribution of large-scale 
material is unconstrained, and thus we do not calculate ex- 
pected ratios for these models. 



